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A multi-block infrastructure for three-dimensional time-dependent numerical relativity 
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We describe a generic infrastructure for time evolution simulations in numerical relativity using 
multiple grid patches. After a motivation of this approach, we discuss the relative advantages of 
global and patch-local tensor bases. We describe both our multi-patch infrastructure and our time 
evolution scheme, and comment on adaptive time integrators and parallelisation. We also describe 
various patch system topologies that provide spherical outer and/ or multiple inner boundaries. 

We employ penalty inter-patch boundary conditions, and we demonstrate the stability and accuracy 
of our three-dimensional implementation. We solve both a scalar wave equation on a stationary 
rotating black hole background and the full Einstein equations. For the scalar wave equation, we 
compare the effects of global and patch-local tensor bases, different finite differencing operators, and 
the effect of artificial dissipation onto stability and accuracy. We show that multi-patch systems can 
directly compete with the so-called fixed mesh refinement approach; however, one can also combine 
both. For the Einstein equations, we show that using multiple grid patches with penalty boundary 
conditions leads to a robustly stable system. We also show long-term stable and accurate evolutions 
of a one-dimensional non-linear gauge wave. Finally, we evolve weak gravitational waves in three 
dimensions and extract accurate waveforms, taking advantage of the spherical shape of our grid lines. 
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I. INTRODUCTION 

Many of the spacetimes considered in numerical rel- 
ativity are asymptotically flat. An ideal kind of domain 
for these has its boundaries at infinity, and at the same 
time has to handle singularities which exist or develop 
within the domain. In a realistic setup, the outer bound- 
ary is either placed at infinity, which is topologically a 
sphere, or one introduces an artificial outer boimdary 
at some large distance from the origin. In both cases, a 
sphere is a natural shape for the boundary. 

There are several possible ways to deal with singular- 
ities. One of the most promising is excision, which was 
first used by J. Thornburg fl"], where he acknowledges 
W. G. Unruh for the idea. Excision means introducing 
a inner boundary, so that the singularity is not in the 
computational domain any more. If done properly, all 
characteristic modes on this inner boundary are leaving 
the domain, so that no physical boundary condition is 
required. Seidel and Suen 12] applied this idea for the 
first time in a spherically symmetric time evolution. 

A well posed initial boundary value problem re- 
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quires in general smooth boundaries 01 • Using spher- 
ical boundaries satisfies all conditions above. Spherical 
boundaries have not yet been successfully implemented 
in numerical relativity with Cartesian grids. However, 
there were many attempts to approximate spherical 
boundaries. For example, the Binary Black Hole Grand 
Challenge Alliance used blending outer boundary condi- 
tions [4, 5], where the blending zone was approximately 
a spherical shell on a Cartesian grid. Excision bound- 
aries often approximate a sphere by having a Lego (or 
staircase) shape |6||, |3, l^- Some of the problems en- 
countered with excision are attributed to this staircase 
shape. A spherical boundary would be smooth in spher- 
ical coordinates, but these are undesirable because they 
have a coordinate singularity on the z axis. A multi- 
block scheme allows smooth spherical boundaries with- 
out introducing coordinate singularities. 

Using multiple grid patches is a very natural thing to 
do in general relativity. When one starts out with a man- 
ifold and wants to introduce a coordinate system, then 
it is a priori not clear whether a single coordinate sys- 
tem can cover all the interesting parts of the manifold. 
One usually introduces a set of overlapping maps, each 
covering a part of the manifold. After discretising the 
manifold, one arrives naturally at multiple grid patches. 

Methods using multiple grid patches in numerical rel- 
ativity were pioneered in 1987 by J. Thornburg 0, 
where he also introduced excision as inner boundary 
condition for black holes. Gomez et al. |10] use two 
overlapping stereographic patches to discretise the an- 
gular direction using the eth formalism. This was later 
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used by Gomez et al. ITlll to evolve a single, non- 
stationary black hole in a stable manner in three di- 
mensions with a characteristic formulation. Thornburg 
iT3l evolves a stationary Kerr black hole in three dimen- 
sions using multiple grid patches using a Cauchy for- 
mulation. Kidder, Pfeiffer, and Scheel I13ll have devel- 
oped a multi-patch pseudospectral code to evolve first- 
order hyperbolic systems on conforming (neighboring 
patches share grid points), touching, and overlapping 
patches. Scheel et al. fT3| used this method with mul- 
tiple radial grid patches to evolve a scalar field on a 
Kerr background. Kidder et al. IT^ used this method 
with multiple radial grid patches to evolve a distorted 
Schwarzschild black hole. 

In this paper, we describe a generic infrastructure for 
time evolution simulations in numerical relativity using 
multiple grid patches, and we show example applica- 
tions of this infrastructure. We begin by defining our 
notation and terminology in section HH where we also 
discuss various choices that one has to make when us- 
ing multi -pat ch systems. We describe our infrastructure 
in section|in]and the patch systems we use in section lTVl 
We test our methods with a scalar wave equation on flat 
and curved backgrounds and with the full vacuum Ein- 
stein equations in sectionslVland rVll We close with some 
remarks on future work in section lVm 



II. MULTI-PATCH SYSTEMS 

Our main motivation for multi-patch systems is that 
they provide smooth boundaries without introducing 
coordinate singularities. This allows us to implement 
symmetric hyperbolic systems of equations for well 
posed initial boundary value problems. However, 
multi-patch systems have three other large advantages 
when compared to using a uniform Cartesian grid. 

No unnecessary resolution. In multi-patch systems the 
angular resolution does not necessarily increase with ra- 
dius. An increasing angular resolution is usually not re- 
quired, and multi-patch systems are therefore more effi- 
cient by a factor 0{n'^) when there are 0{n^) grid points. 

No CFL deterioration for co-rotating coordinates. When 
a co-rotating coordinate system is used, the increasing 
angular resolution in Cartesian grids forces a reduction 
of the time step size. Near the outer boundary, the co- 
rotation speed can even be superluminal. This does not 
introduce any fundamental problems, except that the 
time step size has to be reduced to meet the CFL crite- 
rion. This makes Cartesian grids less efficient by a factor 
of 0{n) for co-rotating coordinate systems when there 
are O(n^) grid points. 

Convenient radially adaptive resolution. In multi-patch 
systems, it is possible to vary the radial resolution with- 
out introducing coordinate distortions 1121] . This is not 
possible with Cartesian coordinates. Fish-eye coordi- 
nate transformations |16] lead to distorted coordinate 
systems. In practise, there is a limit to how large a fish- 



eye transformation can be, while there is no such limit 
for a radial rescaling in multi-patch systems. 

These advantages are so large that we think that fixed 
mesh refinement methods may even be superfluous for 
many applications in numerical relativity. Mesh refine- 
ment can be used adaptively, e.g. to resolve shock waves 
in a star. It would be difficult to handle this with adap- 
tive patch systems. However, mesh refinement is also 
used statically to have higher resolution in the centre 
and lower resolution near the outer boundary. This 
case is very elegantly handled by multi-patch systems. 
Multi-patch systems could also be used to track moving 
features, such as orbiting black holes. 



A. Terminology 

Methods using multiple grid patches are not yet 
widely used in numerical relativity, and this leads to 
some confusion in terminology. 

The notions of multiple patches, multiple blocks, or 
multiple maps are all virtually the same. In differential 
geometry, one speaks of maps covering a manifold. In 
computational physics, one often speaks of multi-patch 
systems when the patches overlap, and of multi-block 
systems if they only touch, i.e., if the blocks only have 
their boundaries in common. However, when discre- 
tised, there is an ambiguity as to what part of the do- 
main is covered by a grid with a certain resolution. See 
figureQlfor an illustration. 

In the following, we say that a grid extends from its 
first to its last grid point. This is different from Thorn- 
burg's notation 1 ,121 : He divides the grid into interior and 
ghost points.^ The interior points are evolved in time, 
and a small number of ghost points are defined through 
an inter-patch boundary condition, which is interpola- 
tion in his case. He defines the grid extent as ranging 
from the first to the last interior point, ignoring the ghost 
points for that definition. Thus, when there are n ghost 
points, he calls "touching" what we would call an "over- 
lap of 2n points". When he speaks of overlapping grids 
according to his definition, then there are parts of the do- 
main covered multiple times, and these overlap regions 
do not vanish in the continuum limit. 



B. Coordinates and tensor bases 

Although not strictly necessary, it is nevertheless very 
convenient to have one global coordinate system cover- 
ing the whole domain. This makes it very easy to set up 
initial conditions from known analytic solutions, and it 



^ Just to demonstrate that terminology can be really confusing, we 
note that Thornburg's notion of "ghost points" is different from 
what Cactus IT^ calls "ghost points". 
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(a) Touching patches, also 
known as blocks. 
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(b) Touching patches with 
additional boundary points. 
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(c) Overlapping patches with 
boundary points. 

Figure 1: Schematics of the difference between touching and 
overlapping grids for different inter-grid relations. The black 
grid points are evolved in time, the hollow grid points are de- 
termined via a boundary condition. 



also makes it possible to visualise the result of a simu- 
lation. While a true physics-based visualisation would 
not rely on coordinates, currently used methods always 
display fields with respect to a coordinate system. 

Since our calculations concern quantities that are not 
all scalars, we have to make a choice as to how to rep- 
resent these. One elegant solution is to use a tetrad 
(or triad) formalism, where one represents vectorial and 
tensorial quantities by their projections onto the tetrad 
elements. However, this still leaves open the choice of 
tetrad elements. 

In a multi-patch environment there are two natural 
choices for a coordinate basis. One can either use the 
global or the patch-local coordinate system. Both have 
advantages and disadvantages. 

In a way, using a patch-local coordinate basis is 
the more natural choice. Given that one presumably 
knows how to evolve a system on a single patch, it is 
natural to view a multi-patch system just as a fancy 
outer boundary condition for each patch. In this way, 
one would continue to use patch-local coordinates ev- 
erywhere, while the inter-patch boundary conditions 
involve the necessary coordinate transformations. It 
should be noted that these coordinate transformations 
mix the evolved variables, because (e.g. for a vector v^) 
what is Vi on one patch is generally a linear combination 
of all Vj in the other patches' coordinate system. 

Using a global coordinate basis corresponds to defin- 
ing a global triad that is smooth over the entire domain. 
This simplifies the inter-patch boundary conditions sig- 
nificantly, since there is no coordinate transformation 
necessary. Instead, one then has to modify the time evo- 
lution mechanism on the individual patches: 



Let the letters a,h,c . . . denote abstract indices for 
quantities in the patch-local tensor basis, and letters 
i, i,k. . . abstract indices in the global tensor basis. The 
triad that defines the global tensor basis is given by fig, 
which is a set of one-forms labelled by a global index 
i. A vector field v is denoted v" in the patch-local tensor 
basis and v' := e'^v" in the global tensor basis. Note that 
v' transforms as a scalar with respect to the patch-local 
tensor basis, as a change in the patch-local tensor basis 
does not change v' . 

When one calculates partial derivatives, e.g. through 
finite differences, one obtains these always in the patch- 
local coordinate basis. It is then necessary to transform 
these to the global coordinate basis. Since the evolved 
variables are scalars with respect to the patch-local co- 
ordinate basis, their first derivatives are co-vectors, and 
their transformation behaviour is straightforward: 

a, = (1) 

where dx" /dx' = {dx'/dx")'^ = (^o)"^ is the (inverse) 
Jacobian of the coordinate transformation. This means 
that using a global tensor basis effectively changes the 
system of evolution equations. 

Using a global coordinate basis also has advantages in 
visualisation. Visualisation tools generally expect non- 
scalar quantities to be given in a global coordinate basis. 
Often, one also wants to examine certain components 
of non-scalar quantities, such as e.g. the radial compo- 
nent of the shift vector. When the shift vector is given in 
different coordinate bases on each patch, then the visu- 
alisation tool has to perform a non-trivial calculation. 

It should be noted that there are many quantities 
which have a non-tensorial character, such as e.g. the 
quantities di^ij of the formulation introduced in iT3l 
which we use below; d/^ij are partial derivatives of the 

three-metric. The quantities (p and T' of the BSSN for- 
malism 1 191 have an even more complex transformation 

behaviour. ^ = In ^^det(7,y) is the logarithm of a scalar 

density, and Y' = is a partial derivative of a tensor 

density. 

It is also necessary to define the set of characteris- 
tic variables at an inter-patch boundary in an invariant 
manner. One convenient way to do so is again to refer 
to a global coordinate basis. That is, one transforms the 
elements of the state vector to the global coordinate ba- 
sis, and can then define the characteristic variables in a 
unique way. 

Last but not least, there is one more compelling argu- 
ment for using a global coordinate basis to represent the 
state vector. Since one, presumably, already knows how 
to evolve the system within a single patch, it may be 
unwise to place all the new complications that a multi- 
patch system brings into the inter-patch boundary con- 
dition. By changing the evolved system to use a global 
coordinate basis, one simplifies the inter-patch bound- 
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aries significantly, and furthermore one can implement 
and test both steps separately. 

III. INFRASTRUCTURE 

We base our code on the Cactus framework fl^,'2^ us- 
ing the Carpet infrastructure Ii2ln22i1 . Cactus is a frame- 
work for scientific computing. As a framework, the core 
("flesh") of Cactus itself contains no code that does any- 
thing towards solving a physics problem; it contains 
only code to manage modules ("thorns") and let them 
interact. Cactus comes with a set of core thorns for basic 
tasks in scientific computing, including time integrators 
and a parallel driver for uniform grids. (A driver is re- 
sponsible for memory allocation and load distribution 
on parallel machines.) 

By replacing and adding thorns, we have extended 
Cactus' capabilities for multi-patch simulations. The 
mesh refinement driver Carpet can now provide stor- 
age, inter-processor communication, and I/O for multi- 
patch systems as well, and the multi-patch and mesh re- 
finement infrastructures can be used at the same time. 
The definitions of the patch systems (see below) and the 
particular inter-patch boundary conditions are handled 
by additional thorns. 

A. Infrastructure description 

A computation that involves multiple patches re- 
quires implementing several distinct features. We de- 
cided to split this functionality across multiple modules, 
where each module is implemented as a Cactus thorn. 
These are 

Driver (D): The driver is responsible for memory allo- 
cation and load distribution, for inter-processor 
communication, and for I/O. It also contains the 
basic time stepping mechanism, ensuring that the 
application e.g. updates the state vector on each 
patch in turn. 

Multi-patch system (MP): The patch system selects 
how many patches there are and how they are con- 
nected, i.e., which face is connected to what face of 
which other patch, and whether there is a rotation 
or reflection necessary to make the faces match. 
The patch system also knows where the patches 
are located in the global coordinate space, so that 
it can map between patch-local and global coordi- 
nates. 

Penalty boundaries (P): This thorn applies a penalty 
boundary condition to the right hand side (RHS) 
of the state vector of one face of one patch. We 
have described the details in |23]. It calls other 
routines to convert the state vector and its RHS on 
the patch and on its neighbour to and from their 



characteristic variables; it is thus independent of 
the particular evolution system. 

Finite differences (FD): As a helper module, we imple- 
mented routines to calculate high order finite dif- 
ferences on the patches I23I1 . These operators use 
one-sided differencing near the patch boundaries. 

Additionally, we make use of the following features 
that Cactus provides: 

Time integrator (TI): The time integrator calls user- 
provided routines that evaluate the RHS of the 
state vector, advances the state vector in time, and 
calls boundary condition routines. 

Boundary conditions (BC): The boundary condition in- 
frastructure keeps track of what boundary condi- 
tions should be applied to what faces and to what 
variables. It distinguishes between inter-processor 
boundaries (which are synchronised by the driver), 
physical boundary conditions (where the user ap- 
plies a condition of his/her choice), and sym- 
metry boundary conditions (which are determined 
through a symmetry of the computational domain, 
e.g. a reflection symmetry about the equatorial 
plane). We have extended the notion of symme- 
try boundaries to also include inter-patch bound- 
aries, which are in our case handled by the penalty 
method. 

Together, this allows multi-patch systems to be 
evolved in Cactus within the existing infrastructure. Ex- 
isting codes, which presumably only calculate the right 
hand side (RHS) and apply boundary conditions, can 
make use of this infrastructure after minimal changes, 
and after e.g. adding routines to convert the state vec- 
tor from and to the characteristic modes. Existing codes 
which are not properly modular will need to be restruc- 
tured before they can make use of this miilti-patch in- 
frastructure. 

We use the penalty method to enforce the inter-patch 
boundary conditions. The penalty method for finite dif- 
ferences is described in 1 24, 25, 26], and we describe our 
approach and notation in l23ll . 

In order to treat systems containing second (or higher) 
temporal derivatives with our current infrastructure, 
one needs to rewrite them to a form where only con- 
tain only first temporal derivatives by introducing auxil- 
iary variables. This is always possible. Strongly or sym- 
metric hyperbolic systems containing second (or higher) 
spatial derivatives 127, 28, 29, ,30, 31, 3ZI could in princi- 
ple also be handled without reducing them to first order. 
The definition of the characteristic modes has then to be 
adapted to such a formulation, and may then contain 
derivatives of the evolved variables. 

This multi-patch approach is not limited in any way 
to finite differencing discretisations. It would equally 
be possible to use e.g. spectral methods to discretise the 
individual patches (as was done in Ill4l1 .) It may even 
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make sense to use structured finite element or finite vol- 
ume discretisations on the individual patches, using a 
multi-patch system to describe the overall topology. 

B. Initial condition, boundary conditions, and time 
evolution 

We now describe how the initial conditions are set up 
and how the state vector is evolved in time. This ex- 
plains in some more detail how the different modules in- 
teract, and what parts of the system have to be provided 
by a user of this multi-patch infrastructure. Each step 
is marked in parentheses (e.g. "(MP)") with the module 
that performs that step. 

Setting up the initial conditions proceeds as follows: 

1. Initialise the patch system. (MP) The patch system 
is selected, and the location and orientation of the 
patches is determined depending on run-time pa- 
rameters. If desired, the patch system specification 
is read from a file (see below). 

2. Set up the global coordinates. (MP) For all grid points 
on all patches, the global coordinates are calcu- 
lated. At this time we also calculate and store the 
Jacobian of the coordinate transformation between 
the global and the patch-local coordinate systems. 
If necessary, we also calculate derivatives of the Ja- 
cobian. Derivatives may be required to transform 
non-tensorial quantities when a patch-local tensor 
basis is used. The Jacobian can be calculated ei- 
ther analytically (if the coordinate transformation 
is known analytically) or numerically via finite dif- 
ferences. 

3. Initialise the three-metric. (User code) Even when 
the evolved system does not contain the Einstein 
equations, we decided to use a three-metric. This 
is a convenient way to describe the coordinate 
system, which — even if the spacetime is flat — 
is non-trivial in distorted coordinates. For exam- 
ple, polar-spherical coordinates can be expressed 
using the three-metric 7,j = diag (1, r'^, r'^ sin^ 0), 
and doing so automatically takes care of all geom- 
etry terms. 

This step is performed by the user code, and it is 
not necessary to use an explicit three-metric in or- 
der to use this infrastructure. At this time, we ini- 
tialise the metric at the grid point locations, speci- 
fying its Cartesian components in the global tensor 
basis. 

4. Initialise the state vector. (User code) Here we ini- 
tialise the state vector of the evolution system, 
also again specifying its Cartesian components in 
the global tensor basis. When evolving Einstein's 
equations, this step and the previous are com- 
bined. 



5. Convert to local tensor basis (if applicable) (User code) 
It is the choice of the user code whether the state 
vector should be evolved in the global or in the 
patch-local tensor bases. We have discussed the 
advantages of either approach above. If the evolu- 
tion is to be performed in patch-local coordinates, 
then we transform the three-metric and the state 
vector at this time. Note that this transformation 
does not require interpolation, since we evaluated 
the initial condition already on the grid points of 
the individual patches. 

At this stage, all necessary variables have been set up, 
and the state vector is in the correct tensor basis. 

The time evolution steps occur in the following hier- 
archical manner: 

1. Loop over time steps. (D) The driver performs time 
steps until a termination criterion is met. At each 
time step, the time integrator is called. 

(a) Loop over substeps. (IT) We use explicit time 
integrators, which evaluate the RHS multiple 
times and calculate from these the updated 
state vector. 

i. Evaluate RHS. (User code) The RHS of the 
state vector is calculated, using e.g. the 
finite differencing thorn. 

ii. Apply boundary conditions to RHS. (BC) 
We decided to apply boundary condi- 
tions not to the state vector, but instead 
to its RHS. This is necessary for penalty 
boundaries, but is a valid choice for all 
other boundaries as well. In our simula- 
tions, we do not apply boundary condi- 
tions to the state vector itself, although 
this would be possible. Both the inter- 
patch and the outer boundary conditions 
are applied via the multi-patch infras- 
tructure in a way we explain below. 

iii. Update state vector. (TI) Calculate the 
next — or the final — approximation of 
the state vector for this time step. 

iv. Apply boundary conditions to the state vec- 
tor. (User code) In our case, nothing hap- 
pens here, since we apply the boundary 
conditions to the RHS of the state vector 
instead. 

(b) Analyse simulation state. (D) The driver calls 
various analysis routines, which e.g. evaluate 
the constraints, or calculate the total energy 
of the system, or output quantities to files. 

The multi-patch thorn knows which faces of which 
patches are inter-patch boundaries and which are outer 
boundaries. Inter-processor boundaries are handled by 
the driver and need not be considered here. Symmetry 
boundary conditions (such as e.g. a reflection symme- 
try) are currently not implemented explicitly, but they 
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can be trivially simulated by connecting a patch face to 
itself. This applies the symmetry boundary via penal- 
ties, which is numerically stable, but differs on the level 
of the discretisation error from an explicitly symmetry 
condition. 

The boundary conditions are applied in the following 
way: 

1. Loop over all faces of all patches. (MP) Traverse all 
faces, determining whether this face is connected 
to another patch or not. Apply the boundary con- 
dition for this face, which is in our case always a 
penalty condition. 

(a) Apply a penalty to the RHS. (P) Apply a penalty 
term to the RHS of all state vector elements. 
This requires calculating the characteristic 
variables at the patch faces. 

i. Determine characteristic variables. (User 
code) Calculate the characteristic vari- 
ables from the state vector on the patch to 
which the penalty should be applied. If 
applying an inter-patch boundary condi- 
tion, calculate the characteristic variables 
from the other patches' state vector as 
well. If this is an outer boundary, then 
specify characteristic boundary data. 

ii. Apply penalty. (P) Apply the penalty cor- 
rection to the characteristic variables. 

iii. Convert hack from characteristic variables. 
(User code) Convert the characteristic 
variables back to the RHS of the state vec- 
tor. 

The edges and corners of the patches require some 
care. In our scheme, the edges and corners of the patches 
have their inter-patch condition applied multiple times, 
once for each adjoining patch. In the case of penalty 
boundary conditions, the edge and corner grid points 
are penalized multiple times, and these penalties are 
added up. This happens for each patch which shares 
the corresponding grid points. We have described this 
in more detail in |23]. 

For nonlinear equations, there is an ambiguity rn the 
definitions of the characteristic variables and character- 
istic speeds on the inter-patch boundaries. Since we use 
the penalty method, the state vector may be different 
at the boundary points on both sides of the interface. 
When the metric is evolved rn time, then it is part of the 
state vector, and it may be discontinuous across the in- 
terface. The definition of the characteristics depends on 
the state vector in the nonlinear case. It can thus happen 
that the characteristic speeds are all positive when calcu- 
lated at the boundary point on one side of the interface, 
and all negative when calculated on the other side of the 
interface. One has to pay attention to apply the penalty 
terms in a consistent manner even if this is the case. In 
our scheme (as described above), we always calculate 
the characteristic information using the state vector on 



that side of the interface to which the penalty terms are 
applied. This scheme does not prefer either side of the 
interface, but it is not fully consistent for nonlinear equa- 
tions. 

Instead of using penalty terms to apply boundary 
conditions, one could also apply boundary conditions 
in other ways, e.g. through interpolation from other 
patches, or by specifying Dirichlet or von Neumann 
conditions. Our infrastructure is ready to do so, but we 
have not performed a systematic study of the relative 
advantages of e.g. penalty terms vs. interpolation. 



C. Time integration 

It is common in numerical relativity to use explicit 
time integration methods. These limit the time step size 
to a certain multiple of the smallest grid spacing in the 
simulation domain. In non-uniform coordinate systems, 
one has to determine the smallest grid spacing explic- 
itly, since it is the proper distance between neighbour- 
ing grid points that matters, not the grid spacing in the 
patch-local coordinate systems, and the proper distance 
can vary widely across a patch. If the three-metric, lapse, 
or shift are time-dependent, then the proper distances 
between the grid points also vary with time. 

Furthermore, the maximum ratio between the al- 
lowed time step size and the grid spacing depends 
not only on the system of evolution equations; it de- 
pends also on the spatial discretisation operators that 
are used, on the amount of artificial dissipation, and on 
the strength of the penalty terms. While all this can in 
principle be calculated a priori, it is time consuming to 
do so. 

Instead, we often use adaptive time stepping, for ex- 
ample using the adaptive step size control of the Nu- 
merical Recipes chapter 16.2]. One can specify a 
time integration accuracy that is much higher than the 
spatial accuracy, and thus obtain the convergence order 
according to the spatial discretisation. In practise, one 
would rather specify a time integration accuracy that is 
comparable to the spatial accuracy. Adaptive time step- 
ping would also have the above advantages when used 
on a single Cartesian grid. 

We would like to remark on a certain peculiarity of 
adaptive time stepping. With a fixed time step size, 
instabilities manifest themselves often in such a way 
that certain quantities grow without bound in a finite 
amoimt of simulation time. Numerically, one notices 
that these quantities become infinity or nari (not a num- 
ber) at some point when IEEE semantics 1341 are used 
for floating point operations. With an adaptive step size, 
this often does not happen. Instead, the step size shrinks 
to smaller and smaller values, until either the step size 
is zero up to floating point round-off error, or the time 
integrator artificially enforces a certain, very small min- 
imum step size. In both cases, computing time is used 
without making any progress. This case needs to be 



7 



Scaling test 
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Figure 2: Wall clock time vs. numbers of processors for 100 
time steps of a test problem. We keep the load of each proces- 
sor approximately constant at 125000 and 194500 grid points, 
respectively. Our implementation scales up to at least 128 pro- 
cessors. 

monitored and detected. 



D. Parallelisation 

Our current implementation parallelises a domain by 
splitting and distributing each patch separately onto all 
available processors. This is not optimal, and it would 
be more efficient to split patches only if there are more 
processors than patches, or if the patches have very dif- 
ferent sizes. This is a planned future optimisation. 

We have performed a scaling test on multiple proces- 
sors. We solve a simple test problem on a patch system 
with multiple patches and measure the time it takes to 
take 100 time steps. ^ As we increase the number of pro- 
cessors, we also increase the number of grid points, so 
that the load per processors remains approximately con- 
stant. This is realistic, because one chooses the number 
of processors that one uses for a job typically depending 
on the problem size. Figure |3] shows the results of the 
scaling tests for two such problem sizes. We find that 
our implementation scales well up to at least 128 pro- 
cessors, and would probably continue to scale to larger 
numbers. See |35] for a comparison of other benchmarks 
using Cactus and Carpet. 

It would also be possible to distribute the domain 
onto the available processors by giving (at least) one do- 
main to each processor. This would mean that one splits 



^ This test was performed with a 4th order Runge-Kutta integra- 
tor, the scalar wave equation formulated in a patch-local tensor 
basis, a seven-patch system, the D6_5 differencing operators, and 
a Mattsson-Svard-Nordstrom dissipation operator. We varied the 
number of grid points per patch from 65^ to 253^ to keep the load 
per processor approximately constant. See section lVAl below. where 
these details are explained. 



domains when one adds more processors, introducing 
additional inter-patch boundary conditions. Penalty 
inter-patch boundary conditions are potentially more 
efficient than using ghost zones, since they require an 
overlap of only one single grid point. An Jtth order ac- 
curate finite differencing scheme, on the other hand, re- 
quires in general an overlap of 2n grid points. Penalty 
boundary conditions thus require less communication 
between the patches. A disadvantage of this scheme is 
that the exact result of a calculation then depends on the 
number of processors. Of course, these differences are 
only of the order of the discretisation error. 

Such differences are commonly accepted when e.g. el- 
liptic equations are solved. Many efficient algorithms 
for solving elliptic equations apply a domain decompo- 
sition, assigning one domain to each processor, and us- 
ing different methods for solving within a domain and 
for coupling the individual domains. The discretisation 
error in the solution depends on the number of domains. 
For hyperbolic equations that are solved with explicit 
time integrators, it is often customary to not have such 
differences. On one hand, this may not be necessary to 
achieve an efficient implementation, and on the other 
hand, it simplifies verifying the correctness of a parallel 
implementation if the result is independent of the num- 
ber of processors. However, there are no fundamen- 
tal problems in allowing different discretisation errors 
when solved on different numbers of processors, espe- 
cially if this may lead to a more efficient implementa- 
tion. 



IV. PATCH SYSTEMS 

We have implemented a variety of patch systems, 
both for testing and for standard application domains. 
It is also possible to read patch systems from file. 

Simple testing patch systems are important not only 
while developing the infrastructure itself, but also while 
developing applications later on. Since the application 
has to provide certain building blocks, such as e.g. rou- 
tines that convert to and from the characteristic rep- 
resentation, it is very convenient to test these in sim- 
ple situations. Many patches have distorted local co- 
ordinate systems, and it is therefore also convenient to 
have patches with simple (one-dimensional) coordinate 
distortions. This can be used to test the tensor basis 
transformations — keeping in mind that some variables 
will not be tensorial, but will rather be tensor densities, 
or partial derivatives of tensors, with correspondingly 
more involved transformation behaviours. 

We have currently two types of realistic patch systems 
implemented: 

a. Six patches: This system consists of six patches 
that cover a spherical shell, i.e., a region with Tmin < r < 
''max- We use the same patch-local coordinates as in 
and 123]. This system is useful if the origin is not part of 
the domain, e.g. for a single black hole. See figure|3]for 
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Figure 3: A cut in the equatorial plane of six patches, in which 
four patches are visible. The outer and inner domain bound- 
aries are spheres. There is one radial coordinate spanning 
r = const surfaces, and two angular coordinates perpendicular 
to that. The radial coordinate is smooth across patch bound- 
aries. 
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Figure 4: A cut in the equatorial plane of seven patches, in 
which five patches are visible. The outer boundary is a sphere, 
the inner patch is a cube. There is again one radial coordinate, 
but it does not span r = const surfaces and it is not smooth 
across patch boundaries except at the outer boundary. The two 
angular coordinates are the same as in the six-patch system. 



an illustration. 

h. Seven patches: This system consists of one cubic 
patch that covers the region near the origin, and six ad- 
ditional patches that cover the exterior of the cube until 
a certain radius /"max- We use the same patch-local coor- 
dinates as in 1 23], which are derived from the six-patch 
coordinates above. This system is useful if the origin 
should be part of the domain, e.g. for a single neutron 
star. See figure|4]for an illustration. 

In addition to these two types, we have variations 
thereof, e.g. a system consisting of only one of the six 
patches assuming a sixfold reflection symmetry. We 
have individual patch types as generic building blocks, 
and we can glue them together to form arbitrary patch 
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Figure 5: A cut in the equatorial plane of the imported binary 
black hole patch system. The outer and inner boundaries are 
spheres. Near the boundaries, the coordinate system is similar 
to spherical coordinates, i.e., there is one coordinate direction 
perpendicular to and two direction tangential to the boundary. 



systems. 

After seeking input from the computational fluid dy- 
namics community, where multi-block systems are com- 
monly used to obtain body-fitted coordinate systems, 
we decided that setting up patch systems by hand is too 
tedious, and that commercial tools should be used for 
that instead.^ We therefore implemented a patch system 
reader that understands the GridPro data format. 
This is a straightforward ASCII based format which is 
specified in the GridPro documentation, and support for 
other data formats could easily be implemented as well. 

Using GridPro, we could easily import patch systems 
with two holes and 27 patches (for a generic binary black 
hole system; see figure|5j and with e.g. 30 holes and 865 
patches (for demonstration purposes; see figure |6j. An- 
other advantage of a tool like GridPro is that the grid 
points are automatically evenly distributed over the do- 
main, which may be difficult to ensure if the grid is con- 
structed by hand. 

V. TESTS WITH THE SCALAR WAVE EQUATION 

We test our multi-patch infrastructure with a scalar 
wave equation on an arbitrary, time-independent back- 
ground. Since we express the coordinate distortions via 
a generic three-metric, there is — from the point of view 
of the code — no difference between a flat and a curved 
spacetime. A stationary black hole background requires 
non-trivial lapse and shift fimctions, but these are desir- 
able even for flat spacetimes: a non-zero shift makes for 
moving or rotating coordinate systems (implementing 



^ We are indebted to our esteemed colleague F. Muldoon for teaching 
us about the state of the art in grid generation. 
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to the time evolution equations 

dtu = p 



dtp = ^%p + —di 
dtVi = dip 



^[f^^p + u'Wvj 



with 



(5) 
(6) 
(7) 

(8) 



If discretised with operators that satisfy summation by 
parts, this system is numerically stable with respect to 
the energy 



'-2 



V7 



p^ + DL^H'hiVj 



dV 



(9) 
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for |j6| < a. In this case, this system is symmetric hy- 
perbolic, and its characteristic variables and speeds are 
listed in Isil . 

We present below evolutions of the scalar wave equa- 
tion on a flat background with the seven patch system 
and on a Kerr-Schild background with the six patch 
system. We compare the respective benefits of using a 
global or a patch-local tensor basis, and we study the 
behaviour of scalar waves on a Kerr-Schild background 
as a test problem. 



(b) region near the origin 



Figure 6: A cut in the equatorial plane of the imported demon- 
stration patch system with 30 spherical holes in arbitrary po- 
sitions close to the centre. The resolution in the centre is much 
higher than near the outer boimdary, demonstrating how to 
achieve the same effect as fixed mesh refinement with multi- 
ple patches. 



fictitious forces), and a non-unity lapse could be used 
to advance different parts of the domain with different 
speeds in time, which can improve time integration effi- 
ciency (although we did not use it for that purpose). 

We use the notation a for the lapse, j3' for the shift 
vector, jij for the three-metric, j'^ for its inverse, and 
7 = det(7y ) for its determinant. 

We evolve the scalar wave equation 



□m = 

by introducing the auxiliary variables 

p = dtu 
Vi = Dill. 



(2) 



(3) 
(4) 



A. Comparing global and patch-local tensor bases 

We present here time evolutions of the scalar wave 
equation on a flat backgroimd with the seven-patch sys- 
tem. We compare two formulations, one based on a 
global, the other based on a patch-local tensor basis. 
We also apply a certain amount of artificial dissipation 
to the system, which is necessary because our formula- 
tion has non-constant coefficients. We use two differ- 
ent kinds of artificial dissipation, which were introduce 
by Kreiss and Oliger Is^ and by Mattsson, Svard, and 
Nordstrom p^, respectively. 

We set the initial condition from an analytic solution 
of the wave equation, namely a traveling plane wave. 
We also impose the analytic solution as penalty bound- 
ary condition on the outer boundaries. This is in the 
continuum limit equivalent to imposing no boundary 
condition onto the outgoing characteristics and impos- 
ing the analytic solution as Dirichlet condition onto the 
incoming characteristics. The patch system has the outer 
boundary at r = 3, and the inner, cubic patch has the ex- 
tent [—1; +1]. We use the penalty strength ^ = at the 
inter-patch and at the outer boundaries. See fz^] for our 
notation for the penalty terms. 

The traveling plane wave is described by 



This renders the system into a first order form, leading 



u{t,Xi) — A cos 



2n Ikix' +Lot 



(10) 
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with co^ = 5%ikj. We set A = 1 and fc,- = [0.2, 0.2, 0.2], so 

that the wave length is 5/\/3. We construct the solutions 
for p and u,- from u via their definitions JSj and Jl}, and 
evaluate these at t = to obtain the initial condition. 
The flat background has a = 1, /3' =0, and = (5,^ 
in the global tensor basis; the patch-local metric is con- 
structed from that via a coordinate transformation. 

We use dissipation operators that are compatible with 
summation by parts (SBP) finite difference operators. 
Introduced by Mattsson, Svard, and Nordstrom in Js^ . 
we call them "MSN" operators. They are constructed 
according to 



A 



MSN 
2p 



^^h^VY.-'DlB,D„ (11) 



where 2p is the order of the interior derivative opera- 
tor, e is the dissipation strength, h is the grid spacing, 
E is the norm with respect to which the derivative op- 
erator satisfies SBP, Dp is a consistent approximation of 
a pth derivative, and Bp is a diagonal matrix. The scal- 
ing with grid spacing is in contrast to standard Kreiss- 
Oliger (KO) dissipation operators Issll . where the scal- 
ing /z^P"^ is used. Experience has shown that with the 
MSN scaling it is sometimes necessary to increase the 
dissipation strength when resolution is increased in or- 
der to maintain stability. For this reason we have imple- 
mented SBP compatible KO dissipation operators con- 
structed according to 



4KO 



22p+2 



+1 



(12) 



where as before Z, is the norm of the 2pth order accurate 
SBP derivative operator. This yields a dissipation oper- 
ator with KO scaling that has the same accuracy as the 
SBP derivative operator near the boundary (and one or- 
der higher in the interior). The price is having to use a 
slightly wider stencil. 

We use 21 grid points in the angular and in the radial 
directions. The central patch also has 21 grid points in 
each direction. This is a very coarse resolution. We use 
the D6_5 stencil of |23], which is globally sixth order ac- 
curate, and add compatible artificial dissipation to the 
system of both MSN and KO type as described above. 
We choose a dissipation coefficient e ~ 3.0. We use a 
transition region that is 0.3 times the size of the patch. 
The overall system is then sixth order accurate. 

With a patch-local tensor basis and diagonal norm op- 
erators the system is strictly stable, i.e., the numerical er- 
ror is at any given resolution bounded (up to boundary 
terms) by a constant (see also i23il ). while with a global 
tensor basis, a small amount of artificial dissipation is 
required. However, since we are using restricted full 
norm operators for this test, dissipation is required for 
both the patch-local and patch-global case. 

Figure 13 shows the Loo norm of the solution error vs. 
time up to t = 50. The discretisation using MSN dis- 
sipation is unstable for e = 3.0 when a global tensor 
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Figure 7: Comparing global and patch-local tensor bases, and 
Mattsson-Svard-Nordstrom (MSN) and Kreiss-Oliger (KO) 
dissipation operators. The graphs show the Loo norm of the 
solution error vs. time for a coarse resolution on a seven-patch 
system. Some artificial dissipation is necessary to stabilise the 
system, since it has non-constant coefficients. For this partic- 
ular value of the dissipation strength e, using a global tensor 
basis is unstable with the MSN dissipation operators, but sta- 
ble with the KO operators. With higher values of e, the system 
is stable for both dissipation operators. 



basis is used, but it is stable when a local tensor basis is 
used. Larger values of e also stabilise the global tensor 
basis discretisation. For the KO dissipation, a dissipa- 
tion strength e = 3.0 is sufficient to stabilise both the 
local and global tensor basis formulations. Note that the 
error levels are very similar in all cases, showing that 
the main difference between the patch-local and patch- 
global tensor basis implementations is that more dissi- 
pation is necessary in order to stabilise the system. 

Finally, we show a typical shape for solution errors in 
figure |8| This shows the solution errors in the quantity 
u across the +x block for two different resolutions. The 
simulation started with a Gaussian pulse as initial con- 
dition. The graph shows the errors at the time t = 25 
along the a coordinate line for b = 0, c = 0; this coor- 
dinate line is approximately the (p coordinate line, ex- 
cept that it has a kink at the block interfaces. (See figure 
|3]for an illustration.) b = places the coordinate line 
into the equatorial plane, and c = chooses the cen- 
tre of the block in the radial direction. The coarse block 
had 17 X 17 X 141 grid points with the outer boundary at 
R = 15, and the fine block had twice this resolution. The 
simulation was run with the Dg-s operator, and we ex- 
pect 6th order convergence. We have intentionally cho- 
sen rather coarse resolutions in the angular direction. 
Because the SBP stencils of the Dg-s operator are mod- 
ified on 7 grid points near each boundary, convergence 
at the boundary is not obvious from this graph. How- 
ever, we show in (23il that this system converges indeed 
to 6th order. 
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Figure 8: This figure shows typical solution error shapes across 
a block. It shows the error along the n coordinate line of the 
+x patch of a seven-patch system. This coordinate line corre- 
sponds to a ^ coordinate line in spherical coordinates, but has 
a kink at the patch interface; see figure 01 This figure shows 
two resolutions which differ by a factor of two, and the coarse 
error is scaled according to 6th order convergence. The error 
is largest near the block boundary. Because the SBP stencils 
of the D5_5 operator are modified on 7 grid points near each 
boundary, convergence near the boundary is not obvious from 
this graph. We show 6th order convergence for this case in 
El. 



B. Scalar wave equation on a Kerr-Schild background 

We also present time evolutions of the scalar wave 
equation on a Kerr-Schild (40l section 3.3] background 
with six patches, which we use to excise the singularity. 
We choose the mass M = 1 and the spin a = 0.9 for the 
background. We place the inner boundary at r = 1.4 and 
the outer boundary at r = 201.4. 

We use as initial condition «(0, x') = 0, u, (0, x') = 0, 
and a modulated Gaussian pulse 



p{0,x') = YtmAexp 



{r-Rf 



W2 



(13) 



where we choose the multipole I = 2, m = 2, the ampli- 
tude A = 1, the radius R = 20, and the width W = 1. 
We use conventions such that 



15 

32 TT 



sm 



cos^ (p — sin^ (p 



(14) 



We impose i( = 0, p = 0, and c, = with the penalty 
method as outer boundary conditions. This means that 
this condition is imposed onto the incoming character- 
istic modes. Since the inner boundary is an outflow 
boundary, no boundary condition is imposed there. At 
the outer boundary, p is indistinguishably close to zero 
at f = (much closer than the floating point round- 
off error), so that there is no noticeable discontinuity to 
the initial condition. We use again the penalty strength 
S = Q for both inter-patch and outer boundaries. 



We use the patch-local tensor basis for this example. 
We use 21 grid points in the angular directions and 1001 
grid points in the radial direction. We use here — for no 
particular reason — different discretisation parameters. 
It is our experience that the stability of the system does 
not depend on the particular choice of stencil, as long 
as it satisfies summation by parts l23ll . We use here the 
D3_4 stencil, which is globally fifth order accurate, and 
add compatible MSN artificial dissipation to the system. 
We choose a dissipation coefficient e = 0.2, and we do 
not scale the dissipation with the grid spacing h. The 
overall system is then fifth order accurate. We use a 
fixed time step size Af = 0.05 with a fourth order ac- 
curate Rimge-Kutta integrator. 

Figure |5] shows a snapshot of the simulation at f = 
92.2. At that time, the wave pulse has traveled approx- 
imately half the distance to the outer boundary. The 
inter-patch boimdaries are smooth, although the config- 
uration is not axisymmetric. Figure [TUI shows extracted 
wave forms from this simulation for the I = 2, m = 2 
and for the £ = A, m = 2 modes. The 1 = 2 mode 
is present in the initial condition. The I = \ mode is 
excited through mode-mode coupling. Its amplitude is 
lO'' times smaller than the 1 = 2 mode at this time. The 
t = \ mode has converged at this resolution, i.e., it does 
not change noticeably when the resolution is changed. 

We have also simulated initial data consisting of an 
I = 2, m = Q mode. In this case, both the i = 0, m = 
and the £ = 4, m = modes are excited through mode- 
mode coupling. As expected, the mode-mode coupling 
vanishes when the spin of the background spacetime is 
set to zero. 

We determine the complex quasi normal frequency 
CO = lor + icvj from the extracted wave form of the 
£ = 2, m = 2 mode. We fit the wave seen by an observer 
at radius r = 5 to a function 



/(f) = A sin(a;Rf - (^) exp(a;/f). 



(15) 



This fit is performed for the real and imaginary part of 
the complex frequency as well as for the amplitude A 
and a phase (p. TablelUshows the frequencies we obtain 
from our simulations for the £ = 2, m = 2 and m = —2 
modes and we compare the predictions made by per- 
turbation theory 14211 . A detailed study of quasi- 
normal mode frequencies and excitation coefficients of 
scalar perturbations of Kerr black holes is in preparation 

m. 

For reasons of comparison between a code using the 
global tensor basis and one using the patch-local tensor 
basis, we evolved a similar physical system using both 
of these methods. We now choose a spin of a = 0.5 and 
initial data with an ^ = 2, m = angular dependency. 
The frequencies obtained by both codes, together with 
the predictions from perturbation theory are shown in 
table m We find that the choice of tensor basis has little 
influence on the accuracy of the results. 
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Figure 9: The patch system and the scalar wave configuration at f = 92.2, also enlarging the region near the excision boundary 
inside the horizon. The background is a rotating black hole with a = 0.9, the initial condition is an £ = 2, m = 2 multipole. Note 
the large scale difference between the outer and the inner boundary, which is handled "naturally" and without mesh refinement. 
There are no artifacts visible at the inter-patch boundaries. 
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Figure 10: The £ = 2, m = 2 and the £ = i, m = 2 modes, extracted at r = 40. The £ = 4:, m = 2 mode is excited through 
mode-mode coupling. Its amplitude is 10^ times smaller than the £ = 2, m = 2 mode. All other modes with ( < 4 are zero up to 
floating point round-off error. 



Mode 


Spin 


'^peturbation 


'^numerical 


/ = 2 


m = 


2 


0.9 


0.781638 - 0.0692893/ 


0.796527(030) 


- 0.0680891(010)/ 




/ = 2 


m = 


-2 


0.9 


0.387710 - 0.0935902/ 


0.387678(001) 


- 0.0934718(100)/ 




/ = 2 


m = 





0.5 


0.491962 - 0.094630/ 


0.491824(100) 


- 0.0946523(200)/ 


(local) 


1 = 2 


m = 





0.5 


0.491962 - 0.094630/ 


0.492432(001) 


- 0.0944723(300)/ 


(global) 



Table I: Comparison of the scalar quasi normal frequencies obtained with the multi block method using a global and local tensor 
basis, and the values predicted by perturbative methods ', \\. \2\. The given error estimates for the numerical values come from 
the uncertainty induced by the fitting procedure. Details about that can be found in ii^l . 
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VI. EVOLVING THE VACUUM EINSTEIN EQUATIONS 

We evolve the vacuum Einstein equations using the 
symmetric hyperboHc formulation introduced in lT3l . 
It includes as variables the three metric gjj, the extrin- 
sic curvature Kjj, the lapse a, and extra variables de- 
noted by dj^ij and A,. When all the constraints are sat- 
isfied, these are related to the three-metric and lapse by 
^kij = ^kSij ^^'^ = di^lna. The formulation admits 
any of the Bona-Masso slicing conditions while still be- 
ing symmetric hyperbolic. The shift has to be specified 
in advance as an arbitrary fimction of the spacetime co- 
ordinates t and x'. In the tests below, we use a time har- 
monic slicing condition and a time-independent shift. 
The characteristic modes and speeds are listed in 0]. 

Previous 3D black hole simulations using this for- 
mulation were presented in L44,1 , using a low order 
Cartesian code and cubic excision, hi 1431 . constraint- 
preserving boundary conditions for this formulation 
were constructed; this paper then studies the well- 
posedness of the resulting initial-boundary value prob- 
lem, and tests a numerical implementation of those 
boundary conditions in fully non-linear 3D scenarios as 
well, again with a Cartesian code. 

After some initial (and quite lengthy) experiments 
with a patch-local tensor basis, we decided to use a 
global tensor basis instead.^ We find that patch-local 
tensor bases increase the complexity of the inter-patch 
boundary conditions very much, because the character- 
istic decomposition of the field variables needs to be 
combined with the tensor basis transformation at the 
patch boundaries. Converting the partial derivatives 
into the global tensor basis is trivial in comparison. In 
addition to that, analysing and visualising the output 
of a simulation is also made much easier when a global 
tensor basis is used. 



A. Robust stability test 

The robust stability test in numerical relativity con- 
sists of evolving featureless initial conditions to which 
random noise has been added. This test was initially 
suggested in IS li^ and later refined in the so-called 
Mexico tests |I48|]; see also (4^. The first stage of this 
test has a domain that is periodic in all directions, i.e., 
has a topology. The most difficult stage of the test 
has a spherical outer boundary through which noise is 
injected. 

We implement the topology with a single patch, 
using penalty inter-patch boundary conditions to give 
the system a toroidal topology. We also use a six-patch 
topology and set the incoming modes to Minkowski on 



* We are grateful to O. Sarbach for suggesting and insisting on this. 
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Figure 11: Robust stability test for the Einstein equations. This 
graph compares three different resolution of a single patch and 
the six-patch system. Random noise is initially added to all 
variables. At late times, the Hamiltonian constraint does not 
grow with time; the system is strictly stable. The two rims 
marked n = 81 and 6p n = 21 have the highest resolutions; we 
aborted them before they reached t = 2000, which corresponds 
to 1000 crossing times for the single patch. 

both the inner and outer boundary. The single patches 
have outer boimdaries at x' e [— 1;+1], the six-patch 
system has the inner boundary at r = 1.9 and the outer 
boundary at r ~ 11.9. We use a Minkowski spacetime as 
background and add random noise with an amplitude 
of 10~^ to all variables. Since the Minkowski spacetime 
has no intrinsic scale, this amplitude should be com- 
pared to our floating point accuracy of approximately 
lO-i*^. Terms that are quadratic in the noise amplitude 
have then the same order of magnitude as floating point 
inaccuracies. 

In the runs shown below, we choose the penalty pa- 
rameter S = 0.5. We discretise the domain with 21^, Al^, 
and 81'' grid points per patch. We use the Dg_4 deriva- 
tive operator and its associated KO dissipation with a 
parameter e = 0.5. We also use an adaptive Runge- 
Kutta time integrator. For comparison, we also show 
results using a six-patch system with the D5_5 operator, 
using the same dissipation parameter. 

FigurelTTlshows the Loo norm of the Hamiltonian con- 
straint as a function of time. The constraints remain es- 
sentially constant. Note that the constraints do not con- 
verge to zero with increasing resolution, since the ran- 
dom data are different for each run and do not have a 
continuxmi limit. 



B. One-dimensional gauge wave 

One of the most difficult of the Mexico tests jH H^] 
is the gauge wave test. This is a one-dimensional non- 
linear gauge wave, i.e., flat space in a non-trivial coor- 
dinate system. This setup lives in a domain, i.e., it 
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has again periodic boundaries, which we implement ei- 
ther as manifestly periodic boundary conditions or as 
penalty inter-patch boundary conditions. We use the 
slightly modified gauge wave which has the line ele- 
ment 



ds^ = -Hdt^ + Hdx^ + dy^ + dz^ (16) 



with 



H = exp 



A sin 



2n{x - f) 



(17) 



We choose the wave length L to be the size of our do- 
main, and set the amplitude A to 0.5. 

Our domain is one-dimensional with x e [— 1;+1]. 
Different from |48], we place grid points onto the bound- 
aries. We use either 41 or 81 grid points. Figure IT3 
compares the shape of the wave form at late times to 
its initial shape. Our system is stable and very accurate 
even after 1000 crossing times when we use manifestly 
periodic boundary conditions. When we impose peri- 
odicity via penalties, the system is less accurate, and 
these inaccuracies lead to a drift which finally leads to 
a negative gn, which is unstable. With 81 grid points 
per wave length, however, our system is both stable 
and very accurate after 1000 crossing times even with 
penalty boundaries. 

Standard lore says that a second order discretisation 
requires one to use at least 20 grid points per wave 
length. This would seem to make it excessive to use 81 
grid points per wave length with our fifth order scheme. 
However, this is not so. According to Kreiss and Oliger 
js^l, using 20 grid points per wave length introduces a 
10% error for each crossing time, making it impossible 
to evolve meaningfully for 1000 crossing times. Achiev- 
ing a 1% error after 1000 crossing time requires about 
2000 grid points for a second order scheme, and approx- 
imately 73 grid points for a fourth order scheme. Given 
these numbers, using 81 grid points for our fifth order 
scheme is appropriate. 



C. Weak gravitational waves 

We now consider perturbations of a flat spacetime, us- 
ing the Regge- Wheeler (RW) perturbation theory L50J to 
construct an exact solution to the linearised constraints 
of Einstein's equations, which we evolve with the fully 
non-linear equations. We linearise about the Minkowski 
spacetime, i.e., we choose a background spacetime with 
the ADM mass M = 0. 

For simplicity we consider an £ = 2, m = odd 
parity perturbation. The resulting metric in the Regge- 
Wheeler gauge and in spherical coordinates is 

ds^ = - dt^ + dr^ + r^{d0^ + sin^e dcj)^) (18) 
-6Sr^ sin^e cosOdrdcp 
- 6 ^(Y + rY) sin^ 6 cos 6 dt dcj) 
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Figure 12: Results for the gauge wave test case, comparing 
differencing operators and resolutions. With manifestly pe- 
riodic boundaries, the system is both stable and highly ac- 
curate after 1000 crossing times. With periodicity imposed 
via penalty boundary conditions, the system is less accurate, 
and lower resolutions are finally unstable after gn becomes 
negative. However, with sufficient resolution and high order 
derivatives, our system is still highly accurate after 1000 cross- 
ing times. 



where a prime denotes derivative with respect to r, and 
where i5 is a parameter that determines the "strength" of 
the perturbation (not to be confused with the penalty 
term, for which we use the same symbol). This met- 
ric satisfies the linearised constraints for any functions 
Y(f = 0,r) and Y{t = 0,r), and satisfies all the 
linearised Einstein equations if T satisfies the Regge- 
Wheeler equation 



(19) 



It is simple to construct a purely outgoing, exact solution 
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to the previous equation: 



Convergence of initial condition 



3dF 
r du 



(20) 



where F(u) is an arbitrary function of w = r — t. The 
above metric is very similar to the one in Isill . except 
that ours is in the Regge- Wheeler gauge. 

However, we use a different coordinate condition for 
our evolutions. We construct from the previous met- 
ric initial conditions in Cartesian coordinates for gjj, K^, 
and dj^ij. We set the initial lapse to one and evolve it 
through the time harmonic slicing condition, and set the 
shift to zero at all times. 

We now want to check at what point non-linear effects 
begin to have an effect on the constraints. That is, we 
want to find out what resolutions are required to see that 
the full, non-linear constraints do not actually converge 
to zero. As initial condition for the function Y and its 
time derivative we choose 



(21) 



n = 0,r) = -2i^Bexp(-i^l^).(22) 

The non-linear constraint violations should have an ap- 
proximate quadratic dependence on the amplitude pa- 
rameter (called A below). Figure^] quantifies this vio- 
lation. It displays the L2 and Loo norms of the non-linear 
Hamiltonian constraint, measured in local coordinates, 
for different families of initial conditions as a function 
of resolution. We use the seven-patch system described 
in section HVl with the outer boundary at r = 6, while 
the inner, cubic patch has an extent of ±1. The initial 
condition parameters are B = —A, a = 0.6, and tq = 3, 
with amplitudes S = 10"^ 10"^, lO'^, 10"''. We use the 
D5_5 derivative with grid points on each patch for 
N = 21, 41, 81, 161. At the highest resolution, the numer- 
ical constraints reach their non-zero continuum values. 
This highest resolution corresponds to 24 grid points 
per (J in the radial direction. This number is compa- 
rable to the coarsest resolutions that we use in the sim- 
ulations presented below. This means that those sim- 
ulations use constraint-violating initial conditions, and 
cannot expect the constraints to decrease with increas- 
ing resolution. 

The non-linear constraint violation should be approx- 
imately a quadratic function of the amplitude of the per- 
turbation 5, at least for small values of S. Figure El 
shows that this is indeed the case. It displays the Hamil- 
tonian constraint violation in the L2 norm for the high- 
est resolution of the previous figure as a function of the 
amplitude S. The measured slope is 2.0002. Figure ITSl 
shows a sample evolution of this odd parity initial con- 
dition family, using the six-patch geometry. 

In order to evaluate the accuracy of our code we now 
choose an initial condition corresponding to an exact so- 
lution of the outgoing type described above. We do so 
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Figure 13: Norms of the Hamiltonian constraint for initial con- 
ditions with different amplitudes as a function of the number 
of grid points on each patch in each direction. For each ini- 
tial amplitude we show both the L2 and the Loo norm. Since 
the initial conditions only satisfy the linearised constraints, we 
clearly see the non-linear constraint violation at the highest 
resolutions. 
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Figure 14: Hamiltonian constraint violation in the L2 norm as 
a function of the amplitude of the perturbation for the highest 
resolution. This shows the expected quadratic dependence on 
the amplitude. 



by choosing 



F{u) 



exp 



(23) 



In these evolutions the parameters that determine the 
initial condition are a = 1, tq = 30, and amplitude S = 
10-3, with inner and outer boundaries at r = 10 and 
r = 60, respectively. 

We extract the wave forms from the numerical results 
by calculating the Regge- Wheeler function at each grid 
point. We then average over one radial shell. There is 
no need for interpolating to a sphere, which simplifies 
the extraction procedure greatly, and probably also im- 
proves its accuracy. FigurelTSlshows the numerically ex- 
tracted Y(> at r = 40, and compared to its exact value Y, 
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Figure 15: Evolution of weak gravitational waves. This shows 
the component Kxx of the extrinsic curvature in the equatorial 
plane at f = 9.06. The gravitational wave packet started as 
a spherical shell approximately in the middle between the in- 
ner and outer boimdaries, and has then split into two packets 
which travel outwards and inwards, respectively. 



for an evolution using the Dg-^ derivative, with the dis- 
sipation parameter e = 0.05. We used two resolutions 
with 16 X 16 X 1001 and 22 x 22 x 1401 grid points on 
each patch. The agreement is excellent. 



VII. CONCLUSIONS 
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We have motivated the use of multi-patch systems in 
general relativity, and have described a generic infras- 
tructure for multi-patch time evolutions. Their main 
advantages are smooth boundaries and constant angu- 
lar resolution, which makes them very efficient for rep- 
resenting systems requiring high resolution in the cen- 
tre and having a radiative zone far away. They may 
even render fixed mesh refinement unnecessary in many 
cases. 

We use the penalty method for inter-patch boxmdary 
conditions. It would equally be possible to use e.g. in- 
terpolation between the patches. A direct comparison 
of these different approaches would be very interesting. 
Our evolution systems are first order symmetric hyper- 
bolic, but second order systems could be used as well. 
We use this infrastructure with high-order finite differ- 
encing operators, but other discretisations such as e.g. 
pseudo-spectral collocation methods can also be used. 
Our infrastructure is based on Cactus and Carpet and 
runs efficiently in parallel. 

We have discussed the relative advantages of using 
global and patch-local tensor bases, and we have com- 
pared the accuracy and stability of both approaches. We 
suggest that using a global basis is substantially more 
convenient, both in the implementation of the code and 
in the post-processing of the generated output. 



Figure 16: The upper panel shows the Regge-Wheeler func- 
tion Y vs. time, comparing the numerical and the exact so- 
lution at r = 40. The agreement is excellent, especially for 
a three-dimensional, fully nonlinear code. The lower panel 
shows the scaled differences to the exact solution for two res- 
olutions with 16 X 16 X 1001 and 22 x 22 x 1401 grid points. 
This demonstrates fifth order convergence, as should be the 
case for the Dg_4 operator. 

We have tested this irtfrastructure with a scalar wave 
equation on a fixed, stationary background, and with a 
symmetric hyperbolic formulation of the Einstein equa- 
tions. We have shown that our multi-patch system 
with penalty boundary conditions is robustly stable and 
can also very accurately reproduce the nonlinear gauge 
wave of the Apples with Apples tests, which has been the 
most difficult of these tests for other codes. 

Finally, we have simulated three-dimensional weak 
gravitational waves in three dimensions, using the same 
nonlinear code, and have accurately extracted the grav- 
itational radiation. The latter is made especially simple 
since the wave extraction spheres are aligned with the 
numerical grid. 

We believe that multi-patch systems, which provide 
smooth boundaries, will be an essential ingredient for 
discretising well-posed initial boundary value prob- 
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lems. 
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